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We develop a framework to discuss stability of epigenetic states as first exit problems in dynamical 
systems with noise. We consider in particular the stability of the lysogenic state of the A prophage, 
which is known to exhibit exceptionally large stability. The formalism defines a quantative measure 
of robustness of inherited states. 



Epigenetics is concerned with inherited states in living 
systems, which are not encoded as genes, but as the (in- 
herited) patterns of expressions of genes. Modulation of 
gene expression, or functional genomics, underlies a wide 
number of biological phenomena, e.g. efficient use of nu- 
trients available to an organism at a particular time. A 
familiar example of inherited patterns of gene expression 
is cell differentiation in multicellular organisms, which, if 
once established, can be propagated for long times. The 
stability of epigenetic states is important, as it simul- 
taneously enables an organism to maintain a favorable 
state, and to keep the ability to change that state in a 
coordinated manner, if external conditions so dictate. 

Some of the simplest examples of two-state systems in 
biology are found among bacteriophages, DNA viruses 
growing on bacterial hosts. When in the state of being 
stably integrated into the genome of the host, known 
as lysogeny, one set of viral genes is expressed. When 
on the other hand the virus is performing other tasks, 
such as directing translation of viral proteins leading to 
lysis, killing the bacterial host, other sets of genes are 
expressed. The classical example is the lysogenic state of 
phage A in Escherichia coli Upon infection of an 

E. coli cell, either A enters a path leading to lysis, or it 
enters lysogeny, which can then be passively replicated 
for very long times. Indeed, the wild-type rate of spon- 
taneous loss of A lysogeny is only about 10~ 5 per cell 
and generation [9, a life-time of the order of five years. 
Moreover, this number is but a consequence of random 
activation of another part of the genetic system, the SOS 
response involving RecA, and the intrinsic loss rate has 
in several independent experiments been found to be less 
than 10 -7 per cell and generation possibly as low 

as 2 • 10~ 9 0. The rate of mutations in the part of the 
lambda genome involved in lysogeny is between 10 and 
10 -7 per generation |||J. Epigenetics is therefore in this 
system more stable than the genome itself. 

One may recall that E. Schrodinger in "What is life?" 
M starts by imagining that the stability of genetic inher- 



itance stems from dynamic equilibria involving macro- 
scopic numbers of molecules. On the basis of the then 
recent experimental data on the dependence of mutation 
rate on radiation, this hypothesis is discarded in favor 
of a molecular model of genetic memory, fore-shadowing 
the DNA-RNA machinery. Epigenetics, in particular 
A lysogeny as presented below, gives an example where 
Schrodinger's first idea is essentially correct. The number 
of molecules needed to achieve stability, over biologically 
relevant time-scales, is however surprisingly small, only 
in the order of hundreds, and thus rather in the meso- 
scopic than in the macroscopic range. Such a number 
however nevertheless gives fluctuations in gene expres- 
sion 0, and accordingly but a finite stability of many 
states. 

A stable state can be likened to a control switch that 
is on. For A the analogy is quite direct |0,p|: lysogeny 
is maintained by protein molecules and A DNA around 
an operator Or, which consists of three binding sites 
OrI, Or 2 and Or3, overlapping with two promotor sites 
Prm and Pr, see Fig. [l| At the binding sites either one 
of two regulatory proteins CI and Cro can bind. These 
proteins are produced from the corresponding genes cl 
and cro, which abut Or, and which are regulated by CI 
and Cro. Hence transcription of cro starts at Pr, which 
partly overlaps OrI and Or2, while transcription of cl 
starts at Prm, which partly overlaps Or2 and Or3. The 
affinity to the two promoters of RNA polymerase, the 
enzyme which catalyses the production of mRNA tran- 
scripts from DNA, depends on the pattern of Cro and 
CI bound to the operator sites. The rates of production 
of the two proteins are therefore functions of the concen- 
trations of the proteins themselves, and balance decay 
and dilution in a stable stationary state with 200-300 CI 
and, on the average, few Cro per bacterium H. This 
is a logical switch, because if CI concentration becomes 
sufficiently low, the increased activation of cro increases 
Cro concentration and decreases cl activation, so that 
lysogeny is ended and lysis follows. 
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FIG. 1. Right operator complex, Or,, consisting of the three operators OrI, Or2 and Or3. cl is transcribed when Or3 is 
free and Or2 is occupied by CI. cro is transcribed when both Or2 and OrI is free. CI dimers bind cooperatively to OrI and 
Or2. 



The simplest mathematical model which embodies 
Fig. [j] is a set of coupled equations for the time rate of 
change of numbers of CI and Cro in a cell : 



Nci = <MiVci,iVcro) 

A> C ro = 0Cro(-^CI,iVcro) 



where the net production rates are 



<Pci 

4>Cro 



fci(Nci,N C ro) 
fcro(N C l,N Cm ) 



Nci/tgi 



(1) 



(2) 



The production terms /ci and /cro are functions of CI 
and Cro concentrations. With no Cro in the system, 
the curve of fci vs. CI concentration has been ex- 
perimentally measured pT| . As reviewed in fj], these 
measurements are consistent with the best available 
data on protein-DNA affinities O-pTf, dimerization 



constants |15|, initiation rates of transcriptions of the 
genes, and the efficiency of translation of the mRNA 
transcripts into protein molecules. The decay constant 
Tci is proportional to the bacterial life-time, since CI 
molecules are not actively degraded in lysogeny, while 
Tcro is about 30 % smaller |^q| . We remark that there 
is considerably more experimental uncertainty in the 
binding of Cro, both to other Cro and to DNA, than 
the binding of CI, see e.g. p7j . As a minimal model of 
the switch, we take tci and Tc ro from data, and deduce 
fd and /cro at non-zero concentrations of both CI and 
Cro with a standard set of assumed values of all binding 
constants, as done in Q. Such a model is conveniently 
visualized by the phase space plot in Fig. |2[ 

If system (Q) is in lysogeny, i.e. in the stable equilib- 
rium S, it will stay there indefinitely The system leaves 
lysogeny when external perturbations push CI concen- 
tration to the left of the separatrix. In vivo, as sketched 
above, the important perturbation is RecA-mediated 
self-cleavage of CI, as a by-product of the SOS DNA 
repair mechanism when host DNA is damaged. The 
functional purpose of the switch, for the virus, is hence 



to sense if the host is in danger, and, if so, jump ship. 

If the numbers of CI and Cro were macroscopically large, 
then (0) would be an entirely accurate description of the 
dynamics. The numbers are however only in the range 
of hundreds. The actual production process is influenced 
by many chance events, such as the time it takes for a 
CI or a Cro in solution to find a free operator site, or the 
time it takes a RNA polymerase molecule to find and 
attach itself to an available promotor. If in a time in- 
terval At the expected number of molecules produced is 
fAt, then the number produced in an actual realization 
has scatter y/ fAt. As a minimal model of the switch 
with finite- N noise, we therefore consider the following 
system of two coupled stochastic differential equations, 
with two independent standard Wiener noise sources 

(dLO?\dLO? o y. 



dNci = 4>cidt 
dN Clo = 4>crodt 



gci dujf 1 
gem du> 



Cro 



(3) 



We assume that there is an equal amount of finite- AT noise 
in decay as in production, and the two noise amplitudes 
are hence 



9Gi = y/fci + Ncifrci 



<?Cro — V /Cro + N Cro /T Cr 



(4) 



The problem of escape from a stable equilibrium point 
like S under a dynamics like (Q) is a first-exit prob- 
lem in the theory of stochastic processes. As such it is 
solved in Wentzel-Freidlin theory jl8],[l9) . A special case 
of Wentzel-Freidlin theory is well-known from chemical 
physics, namely if the noise amplitudes (gci,3Cro) are 
constant and the drift field is a potential field. If so, the 
problem of escape from S is Kramers' classical problem 
of thermal escape from a potential well [BO 21 1. The 



more general Wentzel-Freidlin problem has both similar- 
ities and differences to Kramers' problem, as we will now 
explain. 
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Phase space plot of lambda lysogeny 

Solid line is the optimal exit path 
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FIG. 2. The phase space plot of the dynamical system (|lj), and the optimal exit path in the stochastic dynamical system (|3j). 
The lysogenic state is identified with a stable equilibrium S at (Nci, Nc ro ) ~ (225, I). The basin of attraction of this equilibrium 
is delimited by a separatrix (basin boundary), which passes through the unstable equilibrium point A at (Noi, Ncro) ~ (56, 101). 
The separatrix crosses the Cl-axis at Nci ~ 33. Also indicated is the most probable exit path (full line with arrow) from S to 
A. Insert shows a blow-up around unstable equilibrium A. Note that the most probable exit path goes into A at a different 
angle compared to the trajectories of (|l]) going out of A. Parameter values are as in Q. 



Proceeding hcuristically, we note that the probability 
of a given realization of the noise in time [0,T] is 



Prob({u, t c i, W p°}J) 

Jo r C i " r 



(5) 



exp 



r Cr , 



dt 



where we have introduced the diagonal elements of the 
diffusion matrix, r CI = and r Cro = g^ ID . 

Of all the realizations that move the system from S 
to A, the most probable is the one that minimizes the 
action functional 



ci 



r 



0Cl) 2 (Ncro 



CI 



Per 



dt (6) 



where the initial position is S, the final position A, and 
the minimization is taken over all paths that go from S 
to A in time T. If A min » 1 it can be proved (see |]|) 
that the most probable exit point from the basin of at- 
traction of S is indeed A, and the rate of exit is to leading 
order 



Rate(exit) cx exp (— A 



mux j 



(7) 



The asymptotic formula M) also contains a prefactor of 
dimension one over time ]18| , ^9[ , which in the case of a 
potential field reduces to the prefactor in Kramers' for- 
mula |^0| . In our case the prefactor is of order once per 
bacterial generation. The optimal exit path obeys the ap- 
propriate Euler-Lagrange equations, being the extremal 
of variations of A. Since the Lagrangian in (^) is not 
explicitly time-dependent, the Hamiltonian 



U = r ( r CI PCI + r Cro Per. 



PCI0CI ~r" PCro^Cro (8) 



is conserved along the path. The momenta (pci,PCro) 
are conjugate to the generalized coordinates (Nci,Nc IO ), 
and the Euler-Lagrange equations are equivalent to 

Hamilton's equations N — and p — — We 
note that in this auxiliary classical mechanical system, 
the diffusion constants play the role of space-dependent 
elements of an inverse mass matrix, while the drift field 
is somewhat similar to a magnetic potential. We also 
note that the energy of the optimal exit path, the value 
of TL along that path, must be non-negative, since the 
drift field vanishes at the two end points. On the other 
hand, we have in general dA/dT = —E, where E is the 
energy and T is the transit time. It hence follows that 
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the optimal exit path is a zero-energy path from S to A 
under the Hamiltonian in (@). 

Fig H shows the optimal fluctuation path. In contrast to 
thermal escape from a potential well, where the optimal 
path is always opposite in direction and equal in size to 
the drift field, there is no obvious simple prescription 
to directly compute the path from S to A in Fig. ||. 
The Hamiltonian analogy however suggest the follow- 
ing numerical procedure, using the relaxation method of 
computing solutions to 2-point boundary problems in an 
ODE 

We first find a natural parameter in the system, typically 
one of the binding constants, and vary that to get close to 
the bifurcation where the stable and unstable equilibria 
(S and A) coalesce. The diffusion parameters Tci and 



Tcro are then practically constant in a neighbourhood 
around both points, while the drift field is small. We 
can then compute a path between the two points at high 
energy (equivalently, at a small transit time T), starting 
from a trial solution, which is a straight path at constant 

speed. In other words, N is taken constant along the trial 

path, andp= ^% is given by (T^ 1 ) ■ N — (f). The energy is 

dN 

then lowered incrementally, and the optimal path at each 
energy found by relaxation, using the previous solution 
as the trial solution at the new energy. A zero-energy 
path can thus be found close to the bifurcation, and, 
by changing back the parameter in small steps, again 
using relaxation to find each new path, a zero-energy 
path can be found at the original parameter value. The 
zero-energy motion in the intermediate neighbourhoods 
of the two points always needs to be taken care of by a 
local calculation, as explained in p9j. 



Stable and unstable equilibria 

Changes upon variation of Cro-OR3 binding 
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Change of Wentzel-Freidlin action 

Escape rate per generation is exponential in the action 
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FIG. 3. Systemic changes due to changes in affinity of Cro to operator site Or3. The standard value is —15.5 kcal/mol. 
Stronger binding energies are investigated for use in numerical procedure (see main text), and to explore robustness of lysogeny 
to parameter changes. Other parameters are as in jf|. a) Location of stable equilibrium (S) and unstable equilibrium (A) as 
affinity is varied. The two equilibria are born in a bifuraction at affinity —17.78 kcal/mol and (iVci, iVcro) ~ (176,3.73). At 
increasing value of binding energy (weaker binding), the equilibria move apart, as indicated by the arrows. Insert shows the 
Nci value of the lysogenic state (stable equilibrium), as function of the affinity, in kcal/mol. b) Wentzel-Freidlin action as 
function of affinity. The escape rate from lysogeny is exponential in the action, with a prefactor of the order once per bacterial 
generation. At the standard value (-15.5 kcal/mol), the predicted rate from the model is hence about once in 10 13 generations. 
We note that the escape rate depends very sensitively on parameters. A change of affinity by 1 kcal/mol to —16.5 kcal/mol 
gives an action of about 2.5, and an expected lifetime of the lysogenic state of only about ten generations. 



There is an emerging consensus in molecular biology 
and biological physics that chemical networks in living 
cells have to be robust [p4|-p6|. That is, they have not 
only to work under some conditions, but should work 
under a wide variety of conditions, possibly even under 
change or replacement of parts of the network, see e.g. 
recent mathematical modelling of cell signaling |27]] , and 
of a genetic control network for cell differentiation in 



Drosophila |28| . For the A phage, robustness of lysogeny 
has been experimentally established for several large 
modifications of the Or complex . The present work 
allows us to quantify robustness of epigenetic states. A 
state only exists at all if deterministic equations like (|l|) 
have a stable equilibrium with the corresponding prop- 
erties. This state is stable for long times, even if the 
number of molecules involved is small, if the action A in 
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([|) and (0) is much larger than unity. The state is finally 
robust, if under a typical change of a changeable param- 
eter //, the state still exists and is stable. This means 
that A/i • sjnj must be significantly less than A, where 
the typical change A/i could be the change in binding 
energy upon a single point mutation, of order 1 kcal/mol. 

In Figs. ||a,b) we examine lysogenic stability as func- 
tion of one parameter, the binding of Cro to Or3. If 
we first disregard the noise, we see that a change of 
affinity by 2.25 kcal/mol brings the stable and unstable 
equilibria together, such that the lysogenic state disap- 
pears altogether. We also observe a sensitive dependence 
of the position of the unstable equilibrium, while the 
number of CI in the stable equilibrium (lysogenic state) 
only changes by 30 %. The lysogenic state therefore looks 
qualitatively similar over this range of parameters. These 
are features of the model embodied by equations (0) only. 
If we then bring in our model of the noise, equation (^), 
we see that the action A changes from more than 30 
to less than 3 when affinity changes by 1 kcal/mol, the 
approximate change of binding energy under a single 
point mutation. Such a change hence suffices to desta- 
bilize the switch over biologically relevant time-scales, 
and the model is therefore not robust to such changes, 
in contradiction to fl. This implies the presence of 
some additional mechanism, in order for robustness to 
prevail. We stress that this lack of robustness is an in- 
herent property of the model, true for all variations of 
parameters that have been put forward to quantitatively 
describe these generally accepted mechanisms of control, 
including the recent report that Cro may in fact bind 
cooperatively to OrI Il7j]. 

In conclusion, we have examined the general problem 
of excape from a stable equilibrium in more than one 
dimension, and demonstrated how this determines the 
stability of states of genetic networks. In contrast to 
Kramers' escape from a potential well, the stability of 
inherited states in such networks is not a mathemati- 
cally and computationally trivial problem. The most 
likely exit path does not go along a steepest decent of a 
potential - there is no potential. Instead, such a path 
can be described as a zero-energy trajectory between 
two equilibria in an auxiliary classical mechanical sys- 
tem. Finding it involves similar numerical problems as 
e.g. computing heteroclinic orbits in celestial mechanics. 
The overall lesson of this study is that an examination 
of equilibria and their bifurcations with changing param- 
eter values allow us to quantify both the stability and 
the robustness of particular states of a genetic control 
system. 
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